Classical Langevin Dynamics for Model Hamiltonians 
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We propose a scheme for extending the model Hamiltonian method developed originally for study- 
ing the equilibrium properties of complex perovskite systems to include Langevin dynamics. The 
extension is based on Zwanzig's treatment of nonlinear generalized Langevin's equations. The pa- 
rameters entering the equations of motion are to be determined by mapping from first-principles 
calculations, as in the original model Hamiltonian method. The scheme makes possible, in princi- 
ple, the study of the dynamics and kinetics of structural transformations inaccessible to the original 
model Hamiltonian method. Moreover, we show that the equilibrium properties are governed by an 
effective Hamiltonian which differs from that used in previous work by a term which captures the 
coherent part of the previously ignored dynamical interaction with the omitted degrees of freedom. 
We describe how the additional information required for the Langevin equations can be obtained 
by a minor extension of the previous mapping. 



This paper is dedicated to Professor Josef Devreese on 
the occasion of his 65th birthday and his formal retire- 
ment. May he continue to have many more productive 
years in theoretical physics. 



I. INTRODUCTION 

Our understanding of the properties of relativity sim- 
ple materials has advanced markedly over the last half 
century. A wealth of experimental data exists with which 
theory, both as formal analysis and as first-principles 
computation, is in quantitative agreement. More re- 
cently, attention has turned to such complex materi- 
als as multication oxides which, even at this advanced 
stage of condensed matter physics and materials sci- 
ence, still pose a compelling challenge. Those materials 
possess an immensely rich phenomenology - antiferro- 
magnetism and high-temperature super-conductivity in 
cuprates; colossal magnetoresistance, charge, spin, and 
orbital ordering in manganates; ferroelectricity, antifer- 
roelectricity, ferrodistortion, and relaxors in perovskites 
and related materials; and, more recently, a colossal 
temperature-independent dielectric constant in single 
crystal CaCu3Ti40i2 (CCTO). First-principles studies of 
these complex oxides, however, have been limited to equi- 
librium properties or responses of perfect crystals and 
to relatively simple multilayer structures. More difficult 
equilibrium questions have been addressed for the per- 
ovskite ferroelectrics by constructing a simplified model 
Hamiltonian the parameters of which are obtained from 
first principles computations. The construction of the 
model Hamiltonians for the perovskite ferroelectrics pro- 
ceeds as follows. Those optical branches which are unsta- 
ble in the high symmetry cubic structure are identified. 1 
From them, local modes are constructed which approxi- 
mate the optimally localized lattice Wannier functions 2 
and which capture the dominant anharmonicities in the 



total energy. These local modes form the basis for the 
construction of an "effective" or "model" Hamiltonian 2 
simple enough to use for Monte Carlo simulations of equi- 
librium structures and properties. 3-5 The parameters of 
the model Hamiltonian are few enough in number to ad- 
mit computation by first-principles methods via a map- 
ping to the model Hamiltonian. 

There are many dynamic or kinetic phenomena of great 
interest which cannot be addressed via the model Hamil- 
tonian scheme at its present level of development because 
the model Hamiltonians are too incomplete to use for 
molecular dynamics. Among these are nucleation and 
growth in first-order phase transitions, nonequilibrium 
kinetics associated with phase transitions, domain wall 
movement in ferroelectrics, dielectric dissipation in insu- 
lators and many others. One such problem is presented 
by the phenomenon of pressure amorphization. This has 
been successfully addressed at the conceptional level by 
augmenting a highly simplified model Hamiltonian with 
Langivan dynamics, 6, showing how pressure amorphiza- 
tion can occur by the emergence of metastable displacive 
disorder. In the present paper, Zwanzig's theory of non- 
linear generalized Langevin equations 8 is used to create 
a classical dynamic extension of the model Hamiltonian 
method suitable for addressing such problems beyond the 
present scope of that method. 

In Langevin dynamics, a system evolves under the in- 
fluence of its internal dynamics; of a relaxation term 
or, more generally, of a memory kernel which captures 
the coherent part of its interaction with a thermal bath, 
and of the remaining stochastic forces from the bath. 
In Section 2, we show how to decompose the Born- 
Oppenheimer Hamiltonian of a material into a model sys- 
tem, a bath, and an interaction between them in the form 
with which Zwanzig starts. In Section 3, we recapitulate 
Zwanzig's development to add clarity and to make the 
present paper more nearly selfcontained. In Section 4, 
we transform Zwanzig's abstract equations into explicit 
Langevin equations of motion for the model coordinates 
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and momenta. We then comment on implementation in 
Section 5 and conclude there that exploring the compu- 
tational feasibility of the scheme would be worthwhile. 
We emphasize in Section 6 that the equilibrium proper- 
ties of the system are governed by a model Hamiltonian 
which differs from that previously obtained by a term 
which captures the coherent consequences of the interac- 
tion of the model degrees of freedom with those omitted 
from the model Hamiltonians and comment on general- 
izing the formalism. 

2. THE SYSTEM PLUS BATH DECOMPOSI- 
TION 

We make the Born-Oppenheimer approximation and 
confine ourselves to classical dynamics for the nuclear 
coordinates. The resulting Hamiltonian TL contains the 
nuclear coordinates and momenta. We divide the nu- 
clear coordinates and momenta into two groups, the sys- 
tem variables x which contribute all the anharmonicity 
in TL and the bath variables y the dependence of TL on 
which is harmonic. We are thus supposing that the x 
are the important structural variables in analogy with 
the model Hamiltonian method, but are retaining the re- 
maining dynamic variables y of the material viewed as 
providing a thermal reservoir to the system variables, x 
and y are column victors with coordinate subvectors Q 
and q, respectively, and momentum subvectors P and p, 
respectively: 

.-(?). .-(;) 

TL can now be decomposed into a system part TL' s (x), 
a bath part TL B (y), an d a system-bath interaction 

W SB (ar,2/), 

H = H , s (x) + H' Sb (x,y) + H' B (y) (2) 

The Q correspond to the local modes of the model Hamil- 
tonian method, and the potential energy part of TL' s (x) 
corresponds to the entire model Hamiltonian of the prior 
schemes 3-5 . 

The bath Hamiltonian TL' B (y) is harmonic, 

'H' B (y) = \y T -K-y, (3) 

where y T is the row vector (q T ,p T )- K is a real, sym- 
metric, positive-definite matrix in the y space with the 
submatrices 

K= { K K q K QP ) (4) 

In all cases contemplated here, there are no velocity- 
dependent forces so that K qp — K pq = 0. The system 



- bath interaction is linear in y but nonlinear in x in 
general, 

n' SB (x, y) = -\ {a T (x).y + y T .a{x)) , (5) 

where a(x) is a column vector in the y space which is in 
general a nonlinear function of x. 

To make contact with Zwanzig's development, TL must 
be rewritten. TL' S (x) contains both kinetic, T S (P), and 
potential, VJ,(Q), energies, 

H' s (x) = T s (P) + V s (Q), (6) 

with V s (Q) anharmonic. We now rewrite TL in the form 
used by Zwanzig 

H(x,y)=H s (x)+H B (x,y) (7) 

so that 

TL B (x, y) = \{y~ a(x)f ■ K ■ (y - a(x)) , (8) 

absorbing TL' SB into TL B . This introduces an additional 
potential energy term into TL' S so that 

TLs(X) = T s (P) + Vs(Q), (9) 

V s {Q) = V s {Q)-\a T {x)-K. a {x). (10) 

Thus a(x) becomes the instantaneous equilibrium value 
of y, and the non-negative term subtracted from V s (x) 
in (9b) softens Vs(x). Thus the potential energy of 
the model system differs importantly from the original 
model Hamiltonian and should yield improved equilib- 
rium properties. The Hamiltonian (7) is the starting 
point of Zwanzig's development. 

3. ZWANZIG'S FORMAL LANGEVIN DY- 
NAMICS 

3.1 Kerner Dynamics 

One can write Hamilton's equations of motion for the 
system variables in vector form 

if) 

Kerner's procedure 9 collects these into a single equa- 
tion through the introduction of the antisymmetric ma- 
trix A operating in the x space 




The equation of motion then becomes 

x = A-V x TL, (13) 



2 



y = B ■ V y H, (14) 

where B is the analog of A in the y space. 

3.2 Motion in the y Space 

Introducing the subscript t to specify the time at which 
the dependent variables are evaluated, (14) becomes 



y t = B ■ K ■ [y t - a(x t )]. 
the formal solution to (15) can be written as 
y t - a(x t ) = e tB ' K [y - a(x )] 

- f dt'e tlB K -V x a{xt-v)--xt-t> 

J O 



(15) 



(16) 



after partial integration subsequent to the direct integra- 
tion of (15). 

3.3 Motion in the x-Space 

The equation of motion (13) for x is, more explicity, 

x t = A ■ V X H S - A ■ V x a T {x t ) ■ K ■ [y t - a(x t )}. (17) 

Substituting (16) into (17) yields a closed equation of mo- 
tion for x which will turn out to be the desired Langevin 
equation, 

x t = A ■ V x H s {x t ) 

+ f dt' A ■ V x a(x t ) T ■K-e t ' B - K ^ x a{xt- V )-x t -t 

J O 

- A ■ V x a(x t ) ■ K ■ e tB - K [ Vo - a{x )]. (18) 

The first line on the right side of (18) is a deterministic, 
conservative, nonlinear force. The second line contains 
a state-dependent, that is nonlinear, memory function 
which gives rise to a dissipative force which is nonlocal 
in time. The third line contains a state-dependent re- 
sponse to a stochastic force. In the next section we show 
the latter to be Gaussian random colored noise. 

3.4 The Random Force 

Define the vector F t in y space as a force 

Ft = -K-e tB - K [y -a{x )]. (19) 

The equation (18) for x t is solved for x t given x , the 
initial condition. One could in principle also fix y a as 
an initial condition and solve Eq. (18) deterministically 
by molecular dynamics. One would then, for each value 
of x D chosen, have to compute the complete trajectories 
to establish the statistical dynamics of the system over 
and over for different y . In doing so, any external cou- 
pling to a reservoir forcing the system to relax towards 



equilibrium is ignored. Zwanzig instead makes the as- 
sumption that the bath is in thermodynamic equilibrium 
at t Q given x a . This implies that the initial state y a has 
a probability distribution 



PROB[y \x ] 



e -fm B (ya,x ) 
J dy' e- HB (y'o' x 'o) 



(20) 



so that y a is a Gaussian random variable (GRV). F t , [19], 
is a superposition of GRV's and is therefore a Gausian 
random variable itself with 

(Ft) = 0, (21) 
{F t F v )^k B TK-e {t - t,)B - K . (22) 



3.5 The Langevin Equation 

We can now interpret [16] as a generalized, nonlinear 
Langevin equation, 

x t = A ■ [V x H s (x t )} + I dt'M{t, t') ■ x' t (23) 

J o 

+V x a T (xt)-F t , 

where 

M(tf) = V x a(x t f ■ K ■ e^'^^VMxt,). (24) 

4. EQUATIONS OF MOTION FOR Q AND P 

4.1 Q and P 

As asserted earlier, we exclude velocity-dependent 
forces, which simplifies the formalism. We also absorb 
the nuclear masses into the definition of the momenta so 
that 



Hs = \p-P + V(Q), 



(25) 



H B = t;P ■ P + (<1 - a q (Q)) T ■ K qq ■ (q - a q (Q)), (26) 



a(x).( a f)) 



V x a(x) 



V a 9 (Q) 




If we define 



L(t) =K-e 



tBK 



(27) 



(28) 



(29) 
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then L qq is its only nonzero submatrix so that 



M(t, t') 



(V Q a 9 (Q t )) T • L qq {t - t') ■ V Q a g (Q t .) 







(30) 

All this results in the much simpler, explicit equations 
of motion 

Qt = Pt, (31) 
Pt = -V Q V{Q t ) - I dt'M pp (t,t') ■ P t , 

J O 

-V Q a q (Q t ) T -F q , (32) 

M pp (t, t') = V Q a q {Q t ) T ■ L qq {t - t') ■ V Q a q {Q t ), (33) 

F gt = -L qq (t)-[q -a g (Q )}, (34) 

(F qt F qt ,) = k B TL qq (t-t'). (35) 

The remaining task is to simplify L qq (t — t'). 

4.2 Bath Normal Modes 

We now transform to the normal modes of K in the y 
space. K and B are now comprised entirely of a set of 
2x2 diagonal submatrices K\\ and Baa along the main 
diagonal with 



and correspondingly 



-1 

1 



(36) 



(37) 



Inserting these into L qq (t — t') yields 

Lq\, q \> (t - t') = LJ X COSUJ X {t - t')6\\> , (38) 
greatly simplifying (35) and leading to 

M PP (t,t') = J dujg(uj)Lj 2 cosuj(t-t')Wpp(uj;t,t'), (39) 

Wp P (w;t,t') = [J2Ku\-Lo)V Q a qX (Q t )V Q a qX (Q t ,)]/g(Lu), 

(40) 

S (w)=X;^-w) (41) 

A 

From (40), g(w) is the total phonon density of states 
of the bath. 

Still more explicitly, let I be the index of the degrees 
of freedom of S. The equations of motion become 



Qit=Pit, (42) 
Pit - -V Qe V(Q) dt'M' u (t, t')P e t', (43) 

Mu> = J dujg(uj)uj 2 cos(jj{t-t')Wu>(uj;t,t'), (44) 
W u ,(u;t,t') - [J25(uj x -Lj)V Qi a qx (Q)V Q ,a qx (Q)}/g(uj), (45) 

A 

(PqxtF^v) = k B TL>lcosux(t - t')d X X'- (46) 
5. COMMENTS ON IMPLEMENTATION 

To focus our attention, we discuss as an illustrative ex- 
ample the dynamics of a first-order structural phase tran- 
sition under pressure from crystalline phase C to crys- 
talline phase C. We suppose that C is the simpler, more 
symmetric structure. The transition is associated with 
the development of a soft branch of the phonon spectrum 
of C. There is a spinodal in the pressure-temperature 
plane at which C becomes locally unstable. Implemen- 
tation of the Langevin dynamics for studying the phase 
transition proceeds in a series of steps which we sketch 
in the following. 

5.1 Decomposition into System Plus Bath 

Consider pressures such that C is unstable at T=0. 
Calculate the harmonic normal mode spectrum of C. Es- 
tablish which are the unstable modes or branches of C. 
Construct local modes from the unstable modes and iden- 
tify their amplitudes with Q. All other normal modes are 
allocated to the bath, and K and the normal modes q x 
are determined. 

5.2 Mapping from H to Hs + He 

One first assumes simple analytic forms for the Q de- 
pendence of V(Q) and a(Q) motivated by the insights 
underlying the decomposition of the material into sys- 
tem plus both. Practical considerations limit both the 
complexity of the Q dependence of the forms assumed 
for V(Q) and a(Q) and the fineness with which the Q 
and q spaces are sampled. V(Q) is then obtained from 
first-principles electronic-structure computations of the 
total energies for the fixed sample of Q values which are 
fully relaxed with regard to the sampled q variables. The 
a(Q) are then simply equal to the relaxed values of the qA, 
according to Eqs. (2)-(9). The parameters of the simple 
analytic forms for V(Q) and a(Q) are then determined by 
a best fit to the numerical data. Thus the computational 
cost of generating the necessary information on which to 
base the Langevin simulations is increased beyond that 
of the previous model-Hamiltonian constructions only by 
the relaxation. 
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5.3 Constructing the Langevin Equation 

Knowledge of V(Q) yields the internal forces 
VqiV(Q). Knowledge of a(Q) yields the state-dependent 
quantities X7Q e a\(Q). Knowledge of K yields the fre- 
quency spectrum g(uj), the random force autocorrela- 
tion function, and, together with V Q e a\(Q), Ww(t — 
t'). Thus everything entering the Langevin equations is 
known once the parameters of Hs + H-b are established 
by the mapping of Section 5.2. 

5.4 Solving The Langevin Equation 

Apart from the issue of the number of variables to keep 
track of, there are three essential difficulties in standing 
in the way of implementing numerically the solution of 
the Langevin equations: 1) Each of the two terms in the 
equation for Pi is nonlinear. 2) There is a memory, and 
the noise is colored, which require remembering Q, P, 
and F q \ over an interval of time longer than a relevant 
bath phonon period. 3) There are two quite different 
time scales in the problem, the bath phonon period and 
the characteristic relaxation times implicit in the memory 
function. Issues 1.) and 3.) are encountered in molecular 
dynamics. Issue 1.) is not serious, and issue 3.) has been 
confronted with varying degrees of success in MD. Thus 
the main difficulty is the increase in memory requirement 
associated with issue 2.). Assuming that issue 2.) can 
be overcome, one would start the simulation of the first 
order transition from C to C by drawing an initial value 
of x from the truncated canonical ensemble 

e -/m s (x) 

P{X) = fd X >e-?W ( 4? ) 

with x and x' limited to the basin of attraction of C in 
Vs(Q). One would then choose a time step which permit- 
ted adequate sampling of the memory function and the 
random force autocorrelation function. Forward time in- 
tegration could proceed by a suitable generalization of 
the Verlet algorithm adapted for different-integral equa- 
tions. Updating the random force at each time step need 
not be onerous. In molecular dynamics only P and Q 
need be remembered at each time step, and only \7qV 
needs to be updated. In the current Langevin dynam- 
ics, P, Q W qV,W Qa(Q), and F q all need to be updated. 
Thus, the upper bound to the number of degrees of free- 
dom which can feasibly be incorporated into TLs will turn 
but to be much less than can be incorporated into a 
standard molecular dynamics computation with classi- 
cal pairwise forces. The latter now exceeds 10 6 . Losing 
as many as three to four orders of magnitude in feasi- 
ble cell size would still allow quite complex problems to 
be tackled. We therefore propose that investigating the 
computational complexity of the proposed Langevin dy- 
namics scheme with the goal of establishing its feasibility 



and, if feasible, its size limitations would be highly worth- 
while. 

6. CONCLUSIONS 

The current procedure for generating the "model" or 
"effective" Hamiltonian from first-principles 
calculations 3 " 5 , generates only the potential V^(Q) in 
H' s (x), Eq.(6). By transforming the full Hamiltonian of 
Eq.(2) into the Zwanzig form Eq.(7), we have added an 
additional, softening term to V^(Q) to arrive at a system 
potential Vs(Q), Eq.(lO), which includes all of the coher- 
ent effects of the interaction between the reduced model 
system and the rest of the material. If one constructs 
the master equation for the probability distribution of 
x at t given x Q at t = 0, PROB[x,t — x D ], a general- 
ized Fokker-Planck equation, one finds that PROB[x,t 
— x Q ] asymptotically approaches a unique stationary so- 
lution which is a canonical ensemble for the Hamiltonian 
H s (x), not H' s (x), 

e -/3Hs(x) 

PROB[x,t\ Xo ] - Sdx , e . m . (x .y t^ «>, (48) 

which differs from p(x) in Eq.(47) because there is no 
truncation here. This change is a direct consequence of 
the form of TLb, Eq.(8), which states that the bath vari- 
ables relax towards and oscillate around their instanta- 
neous equilibrium positions a(x). Thus, the use of V^(Q) 
and 7ig(x) for the calculation of equilibrium properties 
is conceptually wrong. 

A quantitative analysis of the errors in the current first- 
principles effective-hamiltonian approach to ferroelectric 
perovskites has recently been carried out by Tinte et al 15 . 
They show that not calculating the fully relaxed Vs(Q) 
and calculating instead the unrelaxed V s (Q) does not, 
for the materials studied, introduce a major error. Nev- 
ertheless, as discussed in Section 5.2, the fully relaxed 
calculation is necessary to generate the interaction ampli- 
tudes a(Q) essential for implementing the Langevin for- 
malism. Tinte et al. report that the major error arises 
from ignoring thermal expansion. The import of that 
conclusion for the present analysis is that anharmonic- 
ity in the bath coordinates q should be included. The 
only practical way to do so without losing the distinc- 
tion between system and bath is through the use of a 
psendo-harmonic approximation 16 to the time evolution 
and statistical distribution of the bath variables y. We 
reserve that generalization of the present formalism for a 
future publication. 

We reiterate the conclusion of Section 5.4 that ex- 
ploring the feasibility of numerical implementation of 
the present Langevin dynamics scheme would be highly 
worthwhile because of the great extension of the reach of 
model Hamilton methods 3-5 ' 10-14 which could result. 

This work was supported by the Center for Piezoelec- 
tric Design under ONR Grant N00014-01-1-0365. The 
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